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SECTION  I 


INTRODUCTION 

Assessment  of  the  potential  impact  of  finite  element  concepts,  applied 
in  the  construction  of  numerical  solution  algorithms  for  computational 
fluid  mechanics,  is  required  and  under  active  study.  The  formal  elegance 
of  the  methodology  has  produced  a  sound  theoretical  basis  for  the  com¬ 
prehensive  computational  simulation  capabilities  now  extant  throughout 
structural  mechanics  (ref.  1).  The  verification  of  this  impact  in  fluid 
mechanics  remains  to  be  achieved,  and  is  the  principal  focus  of  this  re¬ 
search  project. 

In  its  most  elementary  interpretation,  finite  element  theory  returns 
calculus  and  vector  field  theory  to  the  construction  of  discrete  simula¬ 
tion  algorithms  for  any  branch  of  mechanics.  Of  necessity,  using  Taylor 
series  expansions,  one  must  always  be  able  to  verify  the  equivalent  (finite 
difference)  order-of-accuracy,  for  any  of  the  familiar  derivative  terms 
within  the  governing  partial  differential  equation  system.  Specifically, 
linear  (quadratic)  finite  elements  yield,  equivalently,  second- (fourth- ) 
order  finite  difference  representations  for  linear  spatial  derivatives. 

For  other  than  linear  space  derivatives,  however,  the  finite  element  con¬ 
struction  yields  expressions  that  are  not  usually  familiar,  although  upon 
dissection,  can  be  related  to  appropriate  Taylor  series  expansions.  The 
important  feature  is  that  the  theory  produces  the  discrete  analog  expres¬ 
sions,  completely  independent  of  the  a  postereori  ability  to  construct  an 
equivalent  difference  representation. 

In  fluid  mechanics,  confidence  in  the  theoretical  statement  can  only 
be  attained  through  detailed  numerical  assessments  for  progressively  more 
complicated  (and  non-linear)  pertinent  differential  equation  descriptions. 
This  is  the  basic  mission  of  the  University  of  Tennessee  project  in  compu¬ 
tational  fluid  mechanics.  Strict  adherence  to  the  convergence  theory  has 
been  verified  for  linear,  scalar  parabolic  partial  differential  equations 
(ref.  2),  using  linear,  quadratic,  and  cubic  Lagrange,  and  cubic  Hermite 
finite  element  interpolations  within  the  basic  theoretical  statement. 
Importantly,  for  a  non-homogeneous  gradient  boundary  condition,  of  the 
type  omnipresent  in  computational  fluid  mechanics,  this  study  verified 


equal  or  higher-order  convergence  for  the  discrete  analog  produced  by 
the  theory.  Reference  3  documents  adherence  to  convergence  theory,  for 
solutions  of  the  mildly  non-linear  (parabolic)  laminar  boundary  layer 
equations,  for  both  linear  and  quadratic  embodiments  of  the  finite  element 
theory.  Extraneous  error  mechanisms  served  to  obliterate  attainment  of 
the  theoretical  performance  for  the  cubics.  Similar  verifications  for 
the  linear  and  quadratic  formulation,  for  the  consequentially  non-linear 
parabolic  turbulent  boundary  layer  equations  is  reported  in  reference  4. 

Of  primary  theoretical  consideration,  the  Sobolev  norm  used  to  measure 
convergence  was  a  strongly  non-linear  function  of  the  dependent  variable 
set,  yet  the  theory  accurately  quantized  algorithm  performance, 

Viscous  flows  at  large  Reynolds  number,  and  inviscid  flows  are  charac¬ 
terized  by  dominance  of  the  substantial  derivative.  The  generally  dis¬ 
persive  character  of  the  discrete  analog  is  the  dominant  error  mechanism. 
Reference  5  documents  accuracy  and  convergence  assessments  for  scalar 
convection  problems,  driven  by  linear  and  non-linear  hyperbolic  equations. 
The  numerical  results  were  highly  encouraging,  specifically  with  respect 
to  phase  accuracy  and  the  selectivity  of  the  derived  dissipative  finite 
element  algorithm  statement.  The  results  reported  herein  are  an  exten¬ 
sion  of  these  theoretical  and  numerical  considerations.  The  basic  von 
Neumann  stability  analysis  has  been  refined  and  extended.  Numerical  re¬ 
sults  are  presented  for  solution  of  the  complete  Euler  equations,  for 
shocked  one-dimensional  flows,  that  firmly  quantized  performance  and 
resolution  of  discontinuous  variables.  The  formulation  is  extended  to 
a  two-dimensional  description  in  generalized  coordinates,  including  de¬ 
tailed  construction  of  the  tensor  matrix  product  form  of  the  algorithm 
Jacobian.  For  low  speed  flow  prediction,  a  differential  constraint  theo¬ 
retical  formulation  for  the  continuity  equation  is  derived  and  evaluated. 
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SECTION  II 


PROBLEM  STATEMENTS 

The  requirement  is  to  assess  the  key  aspects  of  accuracy,  convergence, 
stability  and  efficiency  of  finite  element  numerical  solution  algorithm 
concepts  applied  to  computational  aerodynamics.  These  performance  measures, 
as  quantized  by  analysis  and  numerical  experiments  on  simplified  scalar 
equations  modeling  key  aspects  of  the  governing  Navier-Stokes  equations, 
were  reported  in  reference  5  .  The  current  research  extends  the  developed 
algorithm  concepts  to  differential  equation  systems  governing  certain 
problem  classes  in  aerodynamics.  Each  class  is  selected  to  permit  isolation 
of  a  key  theoretical  and/or  practical  aspect.  This  section  sumnarizes 
the  various  differential  equation  descriptions. 

The  partial  differential  equation  set  governing  transient,  three- 
dimensional  aerodynamic  flows  is  the  familiar  and  very  non-linear  Navier- 
Stokes  system.  Each  equation  system  studied  is  derived  from  the  Navier- 

Stokes  equations.  In  non-dimensional  conservation  form,  using  Cartesian 

tensor  suirmation  notation,  the  equation  system  governing  flow  of  a  com¬ 
pressible,  viscous,  heat-conducting  fluid  is 

+  O) 

J 

3(puj  ~  r 

L(pui}  a”lt“+HT  uj  pui  +p6ij  “  aij  =0  (2) 

J  L 

L<pe>  '  ^  *  V  '  °1jU1  ‘  V  '  0  (3) 
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In  equations  1-3,  p  is  density,  pu^  is  the  momentum  vector,  p  is  pressure, 
and  e  is  mass  specific  total  energy.  The  Stokes  viscous  stress  tensor  a^. 
and  heat  flux  vector  q.,  in  terms  of  specific  internal  energy  e,  are 
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Assuming  a  polytropic  gas,  p  =  (y-l)pe,  the  convenient  form  for  the  equation 
of  state  is 


L(p)  =  P  -  (y-l)[pe  -  ipUjUj]  =  0  (7) 

Finally,  y  is  the  absolute  viscosity,  k  is  the  coefficient  of  heat  conduc¬ 
tivity,  and  6..  is  the  Kronecker  delta. 

*  J 

One  special  form  of  equations  1-7  for  analysis  corresponds  to  supersonic 
one-dimensional  shocked  flow  in  a  duct  of  variable  cross-sectional  area, 

A(x).  The  specific  form  is 


L(p)  =  +  f^(Pu)  +  PUA"  =  0 

(8) 

L(pu)  =  +  l^tupu  +  p]  +  pu2A'  =  0 

(9) 

L(pe)  =  +  f^Cupe  +  up]  +  pueA'  =  0 

(10) 

where  A'  =  d(jinA)/dx.  Equation  7  is  unchanged  by  limiting  j  =  1.  A 
second  form  for  analysis  describes  two-dimensional,  subsonic  laminar 
flow  of  a  viscous  non-heat  conducting  fluid.  This  is  obtained  from  equations 
1-7  by  setting  k  =  0  and  constraining  the  index  range  1  <,  (i,j)  <  2.  It 
is  further  permissible  to  assume  y  is  constant  for  this  case. 

Many  confined  aerodynamic  flows  exhibit  a  predominant  direction 
of  flow  which  permits  a  simplification  to  equations  1-7  yielding  the 
so-called  parabolic  Navier-Stokes  equations.  Assuming  the  steady,  constant 
density  flow  isothermal  and  turbulent,  and  employing  the  Reynolds  velocity 
decomposition  (ref.  6  ),  the  time-averaged  parabolized  form  of  equations 
1-7  is 
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In  equations  11-13,  the  overbar  signifies  the  time-averaged  mean  variable 
the  underbar  signifies  the  index  is  not  eligible  for  summation,  and  the 
flow  is  assumed  essentially  aligned  with  the  xi  coordinate  direction, 
where  1  j<  (i  ,j)  3. 


A  closure  model 
Present  requirements 


for  the  Reynolds 
are  served  using 


stress  tensor  u j 
the  constitutive 


\l  is  required. 

J 

equation  (ref.  7  ) 
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where  1  <  i  <  3.  The  various  correlation  coefficients  Ca  are  defined 

—  —  v 

(ref.  8).  The  variables  k  and  e  are  the  turbulence  kinetic  energy  and 
isotropic  dissipation  function,  respectively. 


k  =  i  ujuf 


i 

.  2v 

3U^ 

'  3 

Sjk 
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They  are  solutions  to  the  corresponding  parabolic  form  of  the  governing 
differential  equations  (ref.  6), 


L(k)  = 


3X. 


-Ujk  ♦  O-yfc,.  ufuj  |  -  v 


3k_ 

3X. 


3u 


1 


+  u*uj 3^ +  e  =  0 


(17) 
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SECTION  III 


NUMERICAL  SOLUTION  ALGORITHM 


1 .  Finite  Element  Formulation 

The  Navier-Stokes  equation  system,  and  the  various  sub-systems  for 
study  have  been  expressed.  Denote  the  dependent  variable  set  for  any 
given  equation  system  as  {q}  with  members  q..  Then,  for  the  one-dimensional 

T  1 

Euler  equations  7-10,  (q }  *  (p,  pu,  pe,  p},  for  the  two-dimensional 

laminar  flow  Navier-Stokes  simplification  of  equations  1-7,  {q}"*"  =  {p,  pui 
pua,  p},  and  for  the  three-dimensional  turbulent,  incompressible  parabolic 
Navier-Stokes  equations  11-18,  {q }T  =  {ui,  u2,  u3,  P,  k,  e}.  Upon  noting 
that  the  xi  coordinate  in  equations  11-18  spans  the  domain  of  evolution 
of  the  solution,  i.e.,  Xie[xi,x),  the  various  developed  problem  statement 
differential  equation  systems,  excepting  the  incompressible  continuity 
equation  11,  takes  the  general  form. 


3t 


8XJ 


U  .q  •  +  f .  . 
JM1  1J 


+  fi  "  0 


09) 


In  equation  19,  f, . (q . )  and  f . (q  . )  are  specified  non-linear  functions 

I  J  J  *  J 

of  their  arguments,  as  determined  by  the  particular  equation  system. 

For  example,  for  equation  12,  f  . .  =  (p  6.  -  -  a- •  +  u'u^) ,  f .  =  0,  while 
for  equation  8,  f . i  =0  and  f.  =  puA  .° 

The  n-dimensional  partial  differential  equation  system  (19)  is  defined 
on  the  Euclidean  space  Rn,  spanned  by  the  x  coordinate  system  with  scalar 
components  x^,  1  <  i  <  n,  and  t  is  a  generalized  initial-value  coordinate. 
The  solution  domain  fi  is  defined  as  the  product  of  Rn  and  t,  for  all 
elements  of  x  belonging  to  Rn  and  all  elements  of  t  belonging  to  the 
open  interval  measured  from  to,  i.e., 

fi  =  Rn  x  t  =  {(x,t):  xeRn  and  tc[to,t)} 


The  boundary  3fi  of  the  solution  domain  is  the  product  of  the  boundary 
9R  of  Rn,  spanned  by  x,  and  t,  i.e.,  3fi  =  aR  x  t.  Thereupon,  a  differential 
constraint  may  be  applied  of  the  form 
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(20) 


l(qi>  -  aiq,  *  «.  5-  1,-ij  *  «.  -0 

V 

In  equation  20,  the  a.  are  specified  coefficients  and  n..  is  the  outwards 
pointing  unit  normal  vector.  Finally,  an  initial  distribution  for 
on  n0  =  Rn  x  T0  is  required. 

q^x.To)  £  q-(x)  (21) 

The  dissipative  finite  element  solution  algorithm  for  equations  19- 
21  is  a  modest  extension  on  the  form  detailed  in  reference  5  .  The  ap¬ 
proximation  q^(x.,T)  to  the  (unknown)  exact  solution  q(x.,x)  to  equations 
19-21,  is  constructed  from  members  of  a  convenient  finite-dimensional 
subspace  of  Ho(a),  the  Hilbert  space  of  all  functions  possessing  square 
integrable  first  derivatives  and  satisfying  the  boundary  condition  20. 
While  extremely  flexible  in  theory,  the  usual  practice  is  to  employ  the 
most  elementary  functions  (polynomials  truncated  at  degree  k)  defined 
on  disjoint  interior  subdomains  n  ,  the  union  of  which  forms  the  discreti¬ 
zation  of  n.  Hence, 

q^x  ,t)  ;  qJtf.T)  h  "  q*(3,t)  (22) 

e=l 

and  the  elemental  approximation  definition  is 

q*(x,t)  =  (Nk(x)}T{QI(T)}e  (23) 

In  equations  22-23,  i(I)  is  a  free  index  denoting  members  of  {q*1},  and 
sub-  or  super-script  e  denotes  pertaining  to  the  e^  finite  element, 
fle  =  r”  x  t.  The  elements  of  the  row  matrix  { Nk (x) are  assumed  poly¬ 
nomials  on  x.,  1  <  j  <  n,  complete  to  degree  k  and  constructed  to  form 

J 

a  cardinal  basis  (ref.  9). 

The  functional  requirement  of  the  numerical  solution  algorithm  is 
to  render  the  error  in  q^  minimum  in  some  norm.  Based  upon  previous 
experience  (ref.  5  ),  this  is  accomplished  within  the  context  of  the 
finite  element  algorithm  by  requiring  the  error  in  equations  19  and  20, 
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h  h 

i.e.,  L(q")  and  « (qv)  be  orthogonal  to  the  space  of  functions  employed 
to  define  qV.  To  control  non-linearly  induced  instabilities,  it  is  further 

■  L. 

required  that  the  error  in  vLfq1:)  also  be  orthogonal  to  { N.  } -  For  the 

I  K 

parabolic  Navier-Stokes  equations  (as  well  as  the  incompressible  Navier- 
Stokes  equations),  it  is  also  required  that  the  continuity  equation  11 
be  applied  as  a  differential  constraint  on  solution  of  equations  12- 
13.  Identifying  the  (Lagrange)  multiplier  set  6^,  these  four  independent 
constraints  are  linearly  combined  to  form  the  theoretical  statement  of 
the  finite  element  solution  algorithm. 

|  n{Nk)L(qJ)dx  +  Bi-J  n{Nk}vL(q*?)d>£  +  S 2/  n{NR }fi, (q^) dx  +  g3{  nV{Nk}L(p5)dx 
R  R  3R  R 

s  {0}  (24) 

Upon  selection  of  k  in  equation  23,  equation  24  represents  a  system 
of  ordinary  differential  equations  on  t,  of  the  form 

[C]{QI } "  +  [ U] {Q I }  +  [FIJ]{QJ}  +  {FI}  =  {0}  (25) 

A  one-to-one  correspondence  of  terms  in  equations  25  and  19  is  inferred, 
with  the  matrices  in  equation  25  augmented  to  contain  the  various  additional 
terms  introduced  through  g.  f  0  in  equation  24.  (Detailed  expansions 
are  presented  in  a  latter  section.)  An  efficient,  accurate  and  versatile 
integration  algorithm  for  equation  25  is  the  trapezoidal  rule,  see  reference 
2  .  Hence. 

{FI}  =  {QI>i+1  -  (QUj  -  ^[{QD]^  +  {QI)j]  =  {0}  (26) 

defines  the  system  of  non-linear  algebraic  equations  for  determination 
of  the  elements  of  {QI (t) } . 

The  Newton  iterative  solution  algorithm  for  equation  26  is 

[J(FI)]|P+1{6QI>P+11  =  -{FI}P+1  (27) 

The  dependent  variable  in  equation  27  is  the  iteration  vector,  related 
to  the  solution  in  the  conventional  manner. 


(28) 


«»$ s  w>5« +  <«»>!$ 

The  form  of  the  Jacobian  is, 

[J(F1)]  slrgjT  (29) 

the  detailed  construction  of  which  is  direct  upon  expansion  of  the  equation 
24  in  terms  of  the  hypermatrix  formulation  (ref.  3  ). 

2.  Generalized  Coordinates 

The  classical  approach  to  geometric  flexibility  using  a  finite  element 
algorithm  has  been  to  utilize  boundary  conforming  discretizations,  as 
obtained  using  isop&r&setric  triangulation  in  two-dimensions,  for  example, 
and  retaining  the  global  coordinate  system  description  within  equation 
19.  The  alternative  orocedure,  which  has  become  universal  with  finite 
difference  methods  in  computational  aerodynamics,  is  to  generate  a  boundary 
conforming  cGOta^Mte  transformation,  and  manipulate  the  solution  statement 
19  onto  the  transformed  coordinate  system.  Only  modest  differences  exist 
between  this  concept  and  the  isoparametric  finite  element  formulation. 

However,  the  transformation  of  the  algorithm  statement  24  to  the  transformed 
(regular)  grid  introduces  a  considerable  flexibility  and  efficiency  previously 
absent  in  the  classical  finite  element  formulation. 

The  basic  requirement  is  to  generate  and  utilize  a  regularizing 
coordinate  transformation,  that  maps  generally  curved-sided  domain  boundaries 
aRn  onto  coordinate  surfaces  of  the  unit  square  (cube).  The  coordinate 
system  rf  spanning  the  transformed  domain  is  (assumed)  orthogonal  and 
regular;  the  coordinate  transformation  is  simply 

x.-x.fnj.)  (30) 

The  particular  procedure  utilized  to  establish  equation  30  is  ancilliary 
to  this  development,  but  includes  methods  utilizing  Poisson  equation 
solutions  (ref.  10)  or  any  other  analytical  or  numerical  procedure  (cf., 
ref.  11  ).  Specifically,  one  candidate  algebraic  procedure  is  application 
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of  the  classical  finite  element  isoparametric  and  bi-quadratic  cardinal 
basis  {N2(n)l  to  interpolation  of  the  nodal  coordinate  distributions 
defined  for  sub-domains  rJJ  of  Rn.  Figure  1  illustrates  the  concept; 
referring  to  equation  23,  the  transformation  statement  is. 


x.  5  {N2(n)}'{XI}e 


where  the  elements  of  {XI}  are  the  values  that  x..  takes  at  the  location 
of  vertex  and  non-vertex  nodes  defining  the  (curved)  boundaries  of  the 
global  subdomain  of  R^.  For  reference,  the  elements  of  the  cardinal 
basis  (Mrii.nz)}  for  R*  are 


(1  -  m)(l  -  n2 ) (-rii  -  n z 
(1  +  m)(l  -  n2)(  T)l  "  H2 

(1  +  m) (1  +  n2) (  m  +  n2 

i  (1  -  m)  (1  +  n2)  (-na  +  n2 
{N2(n,)>  *  j  2(1  -  ni)(l  -  n2) 

1  4  2(1  +  m)(l  -  n!) 

2(1  -  nDO  +  n2) 

2(1  -  m)(l  +  nl) 


The  elements  of  (N2(n)}  for  R^  are  given  in  reference  1  ,  Ch.  8. 

Returning  to  the  issue,  the  differential  requirement  is  transformation 
of  the  divergence  operator  in  equation  19,  i.e. 

03) 

3Xi  3xi  3rij  [  ' 


The  elements  of  the  inverse  Jacobian  J  1  =  [3n-/3x.]  are  evaluated  as. 


_ _ 1 

3x .  "  det[j] 

J 


u  1  j 

[transformed  cofactor  of  J] 


where 


3X. 

[J]  5  7T1 
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Figure  1 


Physical  Domain 


Transformed  Domain 


a)  Two-Dimensional  Domain  R2 


b)  Three-Dimensional  Domain  R3 


:  Biquadratic  Cardinal  Basis  Coordinate  Transformation 
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is  directly  evaluable  using  equation  30  or  31.  Furthermore,  the  differential 
element  for  equation  24  is 


dx  =  det[J]dn 


(36) 


Consider,  for  example,  equation  19  corresponding  to  the  momentum 
equation  2,  and  using  the  Green-Gauss  form  of  the  divergence  theorem,  the 
first  term  in  equation  24  becomes 

JRn<N)L<puJ)dJ  =  J^N)  det[J]dri 


+  <J  pW 

3R 


,,h  h  ,  h.  h 

Jp  i  H  U  ij 


•  n,det[J]dri 

J 


r  m. 

Rn  9nk 


3n , 


ujpui  +  ""‘ij  '  "fj 


det[J]dri 


(37) 


Identify  the  contravariant  components  of  convection  velocity. 


uj  =  det[J] 


3n. 


i3xjj 


(38) 


with  scalar  components  parallel  to  the  coordinate  system.  Note  that  in 
using  an  algebraic  transformation,  equation  31,  det[J]  cancels  in  equation 
34  yielding 


ujj  =  [Cof.J]  u J  (39) 

where  [Cof.J]  is  the  transformed  co-factor  matrix  of  [J].  Otherwise, 
only  the  elements  3\/3xj  are  known  at  node  coordinates,  and  equation 
38  is  evaluated  directly. 

Using  equations  22-23,  38-39,  and  neglecting  the  surface  integral 
for  the  moment,  equation  37  becomes 


J  {N}L(pu^)dx  =>3  J  _  (DET}!{N}{NHN}T{RHOUI};dn 

rM  i  e  R»i  e  e 


-  J  {UBARK>^{N}  t^—  {N}{N}T{RHOUI}  dn 
nn  e  dr\,  e 


-  J  {ETAKI }^{N}  t§-  {N}{N}T{P}  dri 
Rn  e  ank  e 


+  /  n{ETAKJ>3{N}  — —  {NHN}T{SIGIJ}  dn  (40) 
Dn  e  dru  e 


h  n 

For  equation  40,  since  the  (N.  }  for  q.  are  locally  defined  on  R  ,  the 

K  T  G 

limits  for  the  integrals  correspond,  and  Sg  is  the  operator  projecting 
element  contributions  to  the  corresponding  global  matrices,  equation 
25.  Furthermore,  the  determinant  of  J  is  assumed  interpolated  on  the 
element  domain  R^  using  {N^}  and  the  nodal  values.  Similarly,  ) 

is  recast  as  its  interpolate  {ETAKI (and  also  detJ(snk/3x^)) .  The 
e-subscripted  terms  in  equation  40  are  independent  of  and  can  be  extracted 
from  the  integrand,  leaving  only  products. of  the  polynomials  and  derivatives 
in  {N^}.  These  are  directly  evaluable,  independent  of  the  particular 
choice  for  equation  30,  using  numerical  quadrature.  Hence,  the  stan¬ 
dardized  solution  form,  for  the  first  term  in  equation  24,  as  formed 
from  equation  2,  is 


J  n{N}L(pu5)dx 

R 


{DET }^[M3000]  -f  RH0UI }( 


-  {UBARK}g[M30K0]{RH0UI)e 


-  {ETAKI }g[M30K0]{Ple 


+  {ETAKL}g[M30K0]{SIGIL)e 
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In  equation  41,  the  indices  K  and  L  obey  the  tensor  summation  rule,  I 
is  the  free  index  (for  pu*?),  Sg  is  the  assembly  operator,  and  [M30K0] 
is  the  hypermatrix  equivalent  of  a /3n.  (transformed)  contracted  with  corres- 

J 

ponding  element  distributions.  From  the  standpoint  of  the  coordinate 

transformation,  f DET >e  is  the  nodal  distribution  of  det  [0]  on  R^,  while 

{ETAKI}e  and  {ETAKL are  corresponding  nodal  distributions  of  components 

of  J"1  on  Rn. 

e 

Within  this  generalized  coordinate  framework  for  the  finite  element 
algorithm,  the  grid  and  metric  data  required  for  a  numerical  simulation 
are  the  nodal  distributions  of  J-1  =  [Sn^/aXj],  det  [J]  =  det  [ax^/ari^], 
and  the  elemental  partition  measures  a  .  These  parameters  need  only 
be  sufficiently  smooth  such  that  interpolation  on  the  elemental  domain 
Rn  makes  sense.  Using  a  Poisson  equation-generated  coordinate  transformation, 
for  example,  these  parameters  are  globally  smooth  since  they  are  determined 
on  the  global  domain  Rn.  Alternatively,  the  transformation  parameters 
constructed  on  the  union  of  macro-domains Re  is  further  discretized  into 

n  m  it 

the  union  of  finite  element  domains  R  .  and  J  is  smooth  on  Rm. 

e  m 

3.  Tensor  Matrix  Product  Jacobian 

From  ‘•he  standpoint  of  efficiency,  it  is  desirable  to  construct 
a  tensor  matrix  product  form  of  the  Jacobian  defined  in  equations  27  and 
29.  Such  a  construction  is  readily  accomplished,  provided  the  tensor 
product  cardinal  basis  function  set  (Nk(n)}  is  utilized,  spanning  quadri¬ 
lateral  and  hexahedron  domains  on  R2  and  R3  respectively.  In  this  instance, 
the  Jacobian  matrix  [J(FI)],  equation  27,  is  replaced  by  the  tensor  (outer) 
product  construction 

[J  (FI )]  =*  [Ji]  ®[J2]®[J3]  (42) 

Each  component  [Ja]  is  constructed  from  its  definition,  equation  29, 
assuming  interpolation  and  differention  are  one-dimensional.  The  cor¬ 
responding  formal  ism  in  finite  difference  methodology  is  called  approximate 
factorization  cf.,  ref.  12.  Using  equation  42,  the  solution  statement 
27  becomes 

[Ji]®[J2]®[J3]  (QI}j+]  =  -  (FI  }P+1  (43) 
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Define 


[J2]®[J3]{«qi>S*|  H  -{PI 

[J3]{6QI}j+}  =  -{P2}P^]  (44) 

Then,  the  operational  sequence  for  equation  43  is 

[jxhpi}?;;  -  -  {Fi>j+i 
[J2]{P2}J+]  «  -  {PI >5+1 

CJ3]{«QI>5+?  =  -  {P2}j+i  <45) 


Obviously,  other  permutations  of  the  index  structure  for  [Ja]  could 
be  utilized  for  equations  44-45.  The  key  aspect  of  the  tensor  matrix 
product  Jacobian  is  the  replacement  of  the  very  large  (albeit  sparse) 
matrix  [J],  with  a  block-diagonal  structured  matrices  [Ja].  The  primary 
attribute  is  up  to  several  orders  of  magnitude  reduction  in  computer 
core  storage  requirements  for  the  Jacobian,  as  well  as  significantly 
reduced  CPU  to  construct  the  LU  decomposition  and  perform  the  back  sub¬ 
stitution.  It  must  be  emphasized  that  this  procedure  in  no  way  affects 
the  formation  of  (FI),  equation  27,  wherein  lies  the  accuracy  fe*tvdres 
intrinsic  to  the  finite  element  algorithm  statement,  equation  24.  Compromises 
in  the  construction  of  {FI}  will  invariably  produce  inferior  results 
for  equation  19  non-homogeneous  and  a  vector. 

As  an  example  of  construction  of  [Ja],  consider  only  the  initial- 
value  term  equation  37.  The  corresponding  term  in  the  Jacobian,  equation 
29,  using  equations  25,  26  and  41,  is 


9 { FI }  _ 
3{QJ)  "  e 


{DET}g[M3000]6jj 


(46) 
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where  Sjj  is  the  (discrete  index)  Kronecker  delta.  Then,  referring  to 
equation  40, 


[J]e  =  Ae  IDET}g[M3000]  -  J  n  lDET}J{Nk(n)HNk(n)HNk(n)>Tdn  (47) 

Re 

Assume  for  simplicity  the  most  elementary  case,  i.e.,  k  =  1,  n  =  2  and 
x  =  i.e.,  identity  coordinate  transformation.  Equation  47  becomes, 

assuming  M  =  B  for  n  =  2 


A o{0NE)I[B3000]  =  A  JB200]  =  J  (Nx  (x)  }{NX  (x)  }Tdx  (48) 

e  e  e  Ra 

Assuming  the  rectangular  element  domain  R*  described  by  measures  i  and  u>, 
the  evaluation  of  equation  48  yields 


[J]e  =  Ae[B200]  =  siw[B200]  =  ~ 


4  2  12' 
4  2  1 
4  2 
(sym)  4 


(49) 


The  tensor  product  construction  for  this  matrix  involves  the  evaluation 
of  equation  48  constrained  to  one-dimension.  Hence,  denoting  M  e  A  for 
n  =  1, 


and 


[o«]e 


a“[A200]  =  J  a{Mxa)HMxa)}Tdxa 
R- 


1  f2 

6  [1  2J 


[J2]e 


'2  1 
1  2 


(50) 


(51) 


assuming  a*  =  i  and  a*  =  w.  By  keeping  track  of  entry  locations  in  [J], 
it  is  easy  to  show  that 


[JJe®  [J2]fi  =  [J]e  (52) 

The  extension  to  more  complicated  terms  in  the  Jacobian,  and  to  k  >  1 , 
n  >  2,  builds  upon  the  elementary  concept. 
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SECTION  IV 
THEORETICAL  ANALYSIS 

In  reference  5  ,  an  elementary  truncation  error  analysis  was  doc¬ 
umented  for  the  one-dimensional  linear  form  of  equation  24,  assuming 
63  =  0.  The  results  of  numerical  experiments  reported  therein  indi¬ 
cated  that  superior  accuracy  accrued  to  definition  of  a  family  of  co¬ 
efficients  6i,  with  distinct  values  utilized  in  each  term  in  equation 
24,  as  produced  for  the  study  equation 


L(u)  =  +  uo  =  0 


Dx 


(53) 


A  comprehensive  von  Neumann  stability  analysis  has  been  completed 
for  equation  53,  assuning  that  3} =  {vxAg,  v2Ae)i,  where  vj  and  v?  are 
employed  for  the  first  and  second  terms  in  equation  53,  respectively. 
The  Fourier  solution  for  equation  53  is 


q(x,t)  =  V  exp  [iu>(x-U0t)] 
The  semi-discrete  Fourier  solution  is 

qh(jAx,t)  =  V  exp  [iio(jAx-rt)] 


(54) 


(55) 


Here,  F  =0  +  i 6 ,  where  a  and  6  are  real  numbers,  i  =  /-I  ,  u>  =  2ti/X  is 
the  wave  number  for  wavelength  X,  and  V  is  the  initial  distribution. 
Using  the  linear  finite  element  basis,  k=l  in  equation  23,  the  algo¬ 
rithmic  statement  (equations  24-25)  for  equation  53  can  be  written  in 
the  finite  difference  recursion  form  (ref.  5  ), 


[Ca]e  {Q>c 


A  « 

T 


(l+3vi)Qj_]  +  +  (l-3vi  )Qj+i 


(56) 


Un  r 


-Qj-1 


Qj+1 


U0V2 


-Qj-1 


*  2»j  -  Vl 


(57) 


assuming  a  uniform  discretization  of  Ra  of  measure  Aa.  The  second  form 
in  equation  57  emphasizes  the  role  of  U0v2  as  a  "viscosity."  Substitut¬ 
ing  equation  55  in  56-57,  and  proceeding  through  the  lengthy  algebra 
yields 


1  '  if  +  TT  +  +  VlV2td2  -  t  +  0(d6)] 


1  -  !r  *  to  +  0(dt> 


t  vf 


d2  -  V  +  0(d6) 


v2(-d  +  ^  +  0(d5 )  +  \>1  (d  -  -y  +  0(d5 ) 


1  +^  +  0(d‘) 


+  v 


d2  -  ^  +  0(d&) 


(58) 


where  d  -  wAx  and  0  indicates  order.  Expanding  the  denominator  and  retain¬ 
ing  all  terms  to  order  d6  yields 


a  s  Ur 


1  -  Vi (vi  -  v2 )d2  + 


(  -1  .  ViV2  3 /  V 

780  ~12~  Vl^Vl  “ v?) 


d4  +  0 ( d 6 ) 


(59a) 


6  U 


(vi  -  v2  )d  - 


V,  2 


12 


-  Vj (Vx-  v2  ) 


d3  +  0(d5 ) 


(59b) 


Setting  \>i  =  v2  =  v  1/>T5  yields  the  results  of  reference  5. 


o  *  Ur 


6  -  U0 


1  +  0(d6) 
-d3 


12 


+  0(d5) 


(60) 
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as  obtained  by  requiring  the  third  coefficient  in  equation  59a  to  vanish. 
The  corresponding  level  of  artificial  diffusion  (a)  is  documented  as  ex¬ 
cessive  for  the  linear  and  nonlinear  example  problems  computed  in  refer¬ 
ence  5.  By  not  requiring  vj  and  v2  be  equal,  enforcing  sixth  order  accur¬ 
acy  in  equation  59a  yields 


V2 


vjd2 

V?d2) 


(61) 


Noting  that  d  =  toAx=  2-rr/n,  since  the  approximate  solution  resolves  dis¬ 
crete  wave  lengths  =  nAa,  equation  61  can  be  solved  for  a  range  of 
vi  >  0  and  n  >  2.  Figure  2  is  a  plot  of  this  solution.  Only  for  large  n 
is  the  relationship  linear,  for  vi  sufficiently  large,  and  all  solutions 
converge  at  vi  =  1//T5"  =  v2. 

The  linear  and  nonlinear  two-dimensional  pure  convection  test  cases, 
originally  reported  in  reference  5  ,  were  computationally  reexamined. 
Figure  3  illustrates  a  typical  final  solution  for  the  linear  rotating 
cone  test  case,  with  its  attendant  loss  in  peak  level  and  trailing  dis¬ 
persion  wake.  Table  1  sunmarizes  solution  inaccuracy  on  these  basis,  as 
well  as  loss  of  symmetry,  for  a  range  of  vi  and  v2.  The  best  accuracy 
accrues  to  use  of  vj.  «  0.01  and  v2/vi  ~  0.75.  Both  levels  are  well  be¬ 
low  the  optimal  order-of-accuracy  determination,  see  Figure  2,  or  the 
elementary  analysis  results  v2  =  v1  =  v  =  1//T5.  This  determination 
changes  only  modestly  over  the  Courant  number  range  0.1  <  C  <  0.7,  and  is 
unchanged  using  the  conservative  form  for  equation  53. 

Figure  4  illustrates  a  typical  final  solution  for  the  nonlinear 
traveling  square  wave  test  case,  wherein  U0  in  equation  53  is  replaced  by 
the  dependent  variable  q.  Table  2  sunmarizes  solution  accuracy  in  the 
wave  center  region,  in  terms  of  plateau  level,  depth  of  the  precurser  wake 
and  the  spread  of  the  wave  front  (in  units  of  mesh  measure  A).  For  the 
nonconservative  form  of  equation  53,  the  "optimal"  v2  =  1//T5  and  setting 
vi  zero  yields  an  acceptably  accurate  solution,  except  for  the  spread  to 
4A  of  the  original  1A  interpolation  of  the  wave.  Decreasing  v2  by  J? 
sharpens  the  front  and  induces  a  modest  peak  at  the  plateau  interface. 

For  the  conservative  equation  statement. 
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L°9(V2) 


Figure  3.  Final  Station  Solution  for  Rotating  Cone  Test  Case 
Linear  Finite  Element  Algorithm,  vi  =0  =  v2. 


Figure  4.  Final  Station  Solution  for  Non-Linear  Square  Wave 
Test  Case,  Linear  Finite  Element  Algorithm, 
vi  =  .15//TS,  v2  =  2//TS. 


Table  1 

Rotating  Cone  Solution  Accuracy  Summary 
Linear  Finite  Element  Algorithm,  C=0.5 


Dissipation 

Levels 

.  Solution 

Error 

(Percent) 

Vl 

v2 

Peak 

Symmetry  Wake 

.0 

.0 

-7.0 

-1.0 

BB 

.0 

.0076 

-22.0 

-5.0 

SB 

.012 

.0076 

+2.0 

-2.0 

-21.0 

.012 

.0092 

-2.0 

±0.0 

-19.0 

.012 

.0114 

-7.0 

±0.0 

-16.0 

.015 

.0076 

+10.0 

-6.0 

-26.0 

.030 

.0076 

+17.0 

-50.0 

-74.0 

.258 

.258 

-22.0 

- 

0. 

Table  2 

Non-Linear  Square  Wave  Solution  Accuracy  Summary 
Linear  Finite  Element  Algorithm,  C  =  0.5 


Dissipation 

Solution  Error 

Level  s  x  (1//15) 

Wave 

Precursor 

Wave 

Algorithm 

Form 

ma 

V2 

Peak 

Wake 

S(,|X?d 

.0 

1//2 

+  .07 

-.05 

2 

Non-Conservative 

.0 

1.0 

-.04 

-.05 

4 

H 

.0 

/2 

-.13 

-.06 

4 

II 

.0 

1.0 

+  .21 

-.10 

2 

It 

.0 

/2 

+  .09 

-.05 

2 

Conservative 

.0 

2.0 

-.04 

-.05 

4 

it 

.10 

ft 

-.01 

-.05 

2 

ll 

.15 

ll 

.00 

-.05 

2 

11 

.20 

ti 

+  .01 

-.06 

2 

ll 

.10 

1.0 

+  .11 

-.11 

2 

ll 

.1/2 

n 

+  .15 

-.09 

2 

M 

?3 


L(q)  =  It +  (q2)  =  0 

the  "optimal"  level  of  v2  is  multiplied  by  2  for  compensation.  The  wave 
front  sharpness  is  refined  without  inducing  a  plateau  peak  by  setting 
vi/v2  «  0.1.  These  determinations  remain  generally  valid  over  the  range 
of  useful  Courant  numbers. 

The  square  wave  test  was  repeated  for  the  linear  advection  equation 
53.  Even  though  the  problem  specification  is  linear,  the  resultant  ccarse 
grid  introduces  unacceptable  dispersion  error  to  the  non-dissipative  solu¬ 
tion,  see  Figure  5a.  Figure  5b  illustrates  the  error  control  achieved 
usingv2=  1//T5'  and  -  2/3,  and  Table  3  summarizes  solution  fidelity 

on  the  same  basis  as  the  non-linear  square  wave.  These  results  further 
confirm  the  utility  of  distinct  levels  for  \>i  and  v2. 

Table  3 

Linear  Square  Wave  Solution  Accuracy  Sunmary 
Linear  Finite  Element  Algorithm,  C  =  0.5 


1 

Dissipation 

Solution  Error 

Levels 

x  (1//T5) 

Wave 

Precursor 

Wave 

Vl 

^2 

Peak 

Wake 

(%) 

Spread 

(A) 

.0 

.0 

+  .40 

-3 

2 

.5 

1.0 

.00 

-1 

3 

.67 

1.0 

.00 

-1 

3 

1.0 

1.0 

.26 

-5 

2 
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Viewing  these  results,  it  is  apparent  that  modest  levels  of  vi  cah 
exert  a  profound  impact  on  solution  fidelity  in  both  linear  and  nonlinear 
applications.  The  correction  mechanism  appears  as  inducement  of  a  leading 
phase  error  that  compensates  for  the  lagging  error  intrinsic  to  the  basic 
algorithm.  Selecting  optimal  order-of-accuracy  as  the  criterion  for  de¬ 
termination  of  v>2  does  not  appear  inappropriate  for  the  linear  problem, 
for  which  the  Fourier  analysis  is  exact.  In  concert,  this  analysis 
does  estimate  the  near-optimal  level  of  v2  for  the  nonlinear  problem. 
Viewing  Figure  2,  the  coordinate  v2  =  2//T5,  vi/v2  =0.1,  lies  near  the 
centroid  of  the  distribution  on  n,  indicating  existence  of  the  computa¬ 
tional  compromise. 

The  von  Neumann  stability  analysis  has  been  extended  to  the  quadratic 
element  algorithm  statement,  k  =  2  in  equation  23.  For  a  uniform  dis¬ 
cretization,  with  Ae  =  2Ax,  the  finite  difference  recursion  relation  forms 
for  the  two  terms  in  equation  53,  at  the  vertex  nodes  of  the  discretiza¬ 
tion,  are 


Se[[C]e{Q>e  *  33  [-(1+5vi)Qj_2  +  2 (I+IOvi )Qj+1 


+  2(1-10v!)Qj+1  -0-5vx)qj+J  (63) 


Se  [Uj{Q>e  *  ^  0+2v2)Qj_2  -4(l+4v2)Qj-i  +  28v2Qj 

+  4(l-4v2)Qj+1  -(l-2v2)Qj+2] 

=  ~6~  [^J-2  +  4(^j+l  ‘  Qj+2] 


Qj-2  -  8Qj.!  +  14Qj  -  8Qj+1  +  Qj+2] 


The  second  form  for  equation  64  again  emphasizes  the  role  of  v2U0  as  a 
"viscosity."  Inserting  equation  55  in  63-64,  and  proceeding  through  the 
lengthy  algebra  yields  a  simultaneous  equation  system  for  o  and  6.  Ex¬ 
panding  the  determinant  and  retaining  all  terms  to  0(dc)  yields 


U0  1-4\>iV2  +  (~14  +  184viv2  -  60vj  +  240  viv>2 ) 


T50VlV2  +  3vi 


'viv2  +  -  64  vxv2 


0(d6 ) 


Equation  65  indicates  the  quadratic  element  algorithm  only  second  order 
accurate  for  v>i  =  0  and/or  v2  =  0.  Numerical  experience,  to  be  discussed, 
indicates  =  0  and  v2  *  0(l//T5')  is  preferable  for  accuracy.  The  recur¬ 
sion  relations  for  the  non-vertex  nodes  of  the  quadratic  algorithm  are 
provided  by  equations  56-57,  by  making  the  identities  j  =>  j+1  and  j  =>  j-1 
Hence,  the  interaction  of  these  two  relations,  in  concert  with  vi  =  0  , 
must  reinforce  to  provide  the  excellent  phase  accuracy  associated  with  use 
of  the  quadratic  element  formulation. 

As  will  become  documented,  this  linearized  theoretical  analysis  does 
accurately  estimate  approximate  levels  for  vi  and  v2.  Importantly,  the 
theoretical  statement  of  the  finite  element  algorithm,  equation  24,  is 
confirmed  to  induce  a  phase-selective  dissipation  capable  of  controlling 
non-linearly  induced  instabilities  as  well  as  distributed  phase  celerity. 
The  progression  to  more  complete  equation  systems  is  required  for  further 
quantization. 


V.  DISCUSSION  AND  RESULTS 


1.  One-Dimensional  Compressible  Flow 

Equations  7-10  express  the  conservative  formulation  for  one¬ 
dimensional,  inviscid  compressible  flow  in  a  duct  of  variable  cross 


section.  The  nonconservative  form  is 

L(p)  =  ft  +  “lx  +  Hix  +  puA^  =  0  (66) 

L(u)  =  p§£  +  pu|^  +  ||  +  pu2A'  =  0  (67) 

L(e)  =  +  puff  +  +  ^  +  pue/r  =  0  (68) 

L (p)  =  P  -  (Y-l)p(e-iu2)  =  0  (69) 


For  the  dissipative,  finite  element  algorithm  statement,  equation 
24,  33  =0  is  appropriate  for  either  equation  system.  Denoting  the 
elements  of  {QI }  and  {FI},  1  <  I  <  4,  as  {R,  U,  E,  P}  and  {FR,  FU,  FE, 
FP}  for  clarity,  the  algorithm  statement  equation  (26)  then  becomes,  for 
the  nonconservative  continuity  equation  66,  for  example 


{FR}  =  ([A200]  +  v11[A210]){R}j:+1  [{U}T([A3001  ]  +  [A3100]){R} 

+  Vi{U}T([A3011]  +  A31 1 0] ) {R }  +  {U }T ( [ A40001 ] { A} ) {R}1 


J  j+1  ,j 


=  {0} 


(70) 


The  A-prefix  on  the  matrices  in  equation  70  denote  a  one-dimensional 
evaluation  of  the  corresponding  integrals  in  equation  24.  By  defining 
a  common  denominator,  the  elements  of  each  matrix  are  integers  with 
values  dependent  strictly  on  the  degree  k  of  the  approximation  poly¬ 
nomial,  equation  23.  These  matrices,  for  1  <  k  <  2,  are  given  in 


Appendix  A.  The  Boolean  index  1  indicates  the  location  of  the  spatial 
derivative  within  each  term,  the  0  indicates  simple  interpolation,  and 
vi  are  the  dissipation  parameters  for  the  dependent  variable  p.  The  in¬ 
dex  pair  j+1 , j  denotes  evaluation  of  each  term  in  brackets,  using 
(QI}j+l  and  (QI}j,  and  simming,  while  {  }j+1  j  denotes  {QI}j+1  -  (QI)j, 
see  equation  26. 

The  form  of  equation  70  is  particularly  attractive  for  programming, 
as  well  as  generation  of  a  Jacobian,  and  a  modest  expansion  for  k  =  1 
elements  illustrates  the  notational  structure.  In  particular,  see 
Appendix  A,  for  k  =  1, 

£*00],  5  f  [  1  2  ] 

[A210]e  .  f  [;]  |  ]  (71) 

The  time-derivative  term  in  equation  66,  when  inserted  into  equation  24; 
yields  two  terms 

,R1  <Nk)wdx-  6>jRl  4  <N‘' If dx 

=  S  [f  ,  {Nx }{N1}Tdx  {R}' 

G  Re 

-  v*  Ae(  J-  tNiKNi^dx  {R}' 

1  Jr"  dX  J 

e 

=  SeL  [A200]e{R}'  -  [A210]{R}gJ  (72) 

Inserting  equation  71  into  72,  assuming  Ae  uniform,  and  assembling  (Se) 
the  element  contribution  at  a  common  node  "j"  yields  identically 
equation  56.  If  Ae  was  not  uniform,  then  the  finite  difference  recur¬ 
sion  relation  72  would  have  to  reflect  this,  as  is  already  embedded 
within  the  hypermatrix  form  in  equation  70.  By  the  same  token,  for  Ae 
and  ue  uniform,  equation  57  is  the  finite  difference  recursion  relation 


equivalent  of  the  first  convection  term  in  equation  66, 


f^)  u|^dx 
R 

=  Se[Ae(U}J  [A3001]e{R}e  +  vim]  [A3011]e  {R)e"j 
=  |{U}T[A3001]  +  Vi{U}^  [A301 1 ]j (R }  (73) 

The  remaining  two  terms  in  equation  70  stem  from  and  puA hence 

(fRl‘Vi-S*  ‘>uA')dx  -  S‘JR, 

=  $e[VR}e  [A3001  ]{U}e  +  vJ{R}g  [A301 1  ]{U)e 

+  Ae{U}J [[A40001 ]{A}p] {R}  1 
el  ej 

=  {U}T[A3100]{R}  +  vJ{U)T  [A3110]{R> 

*  {U}T  f[A40001]{A} J{R)  (74) 

(  eJ 

The  two  right  side  forms  for  equation  74  illustrate  transposition  of 
pre-  and  post-contraction  vectors.  The  elements  of  {A}g  are  element 
nodal  values  of  JtnA(x),and  [A40001]  is  of  global  rank  two  and  hyper¬ 
rank  two,  see  Appendix  A. 

The  finite  element  algorithm  statements  for  the  remaining  equations 
67-69  are 

(FU)  =  {R)T([A3000]  +  V2[A3010]){U>j+1 

+  ~  J(R)T[[A40010]{U}  (U)  +  [A201 ]{P) 


+  v| { R }T ( C  A401 10]{U}){U} 

+  { R  }T  ( -l  U}T[  A500001 3  ( A} )  { U  }1  =  {0}  (75) 

-1  j+1  »J 

{FE}  =  {R  }T  ( [  A3000]  +  Vj[A3010]){E}j+1  j 

+  y|wT([A40010]{U}){E}+ {U}T([A3001]  +  [A3100]){P} 

+  V2{R}T([A40110]{U}){E} 

+  {R}T({U}TCA500001]{A}){E} 

J  j+1  »j 

»  {0}  (76) 

{FP}  =  [ A200] { P }  -  (y-1  ){R}T[A3000]{E) 

+  ^R}T([A40000]{U}){U)  =  (0)  (77) 

In  equations  75-76,  the  matrix  [A500001]  is  of  hyperrank  three,  i.e., 
the  elements  are  square  matrices  with  column  matrices  as  elements. 
Equation  77  is  simply  the  interpolation  equivalent  of  an  algebraic 
equation. 

The  next  requirement  is  to  construct  the  Jacobian  of  the  iteration 
algorithm,  equation  27.  From  its  definition,  equation  29,  and  equations 
70,  75-77 

=  [JRR]  =  [A200]  +  v[[A210] 

+  y-  (U}T  [A3001  ]  +  [A3100]  +  v*[A3011] 
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9{FR} 

W 


3{FU} 

W 


<lP 

+  Vi[A3110]  +  [A40001 3{A} 

Jj+1 


[JRU]  =  y-{R}T[[A3lOO]  +  [A3001  ]  +  v*[A3nO] 

+  v^[A3011]  +  [A40001]{A}JM 

[JRE]  =  [0] 


P 

j+1 


[JRP]  =  [0] 


[JUR]  =  1{UP+1  -  Uj }T([A3000]  +  v^[A3010]) 


+  y-(U}T 


[A40001]{U}  +  v*[A41100]{U} 

+  {U}T[A500001]{A}j 


P 

j+1 


[JUU]  =  {R}T(A3000]  +  Vj[A3010]) 


+  ^-{R}1 


([A40010]  +  [A40001 ]  +  v|[A40110] 


+  v2[A40101]){U} 

+  2({U}T[A500001 ]{A}) 
[JUE]  =  [0] 


P 

j+1 


(78) 


[JUP]  =  [A201 ] 


(79) 


[JER]  *|{eP+1  -  Ej}T([A3000]  +  vJ[A3010]) 

+  f-  {E)T[[A40001]{U}  +  v|[A41100]{U} 

+  {U}T[A500001  ]  {A}1 1 


Jj+1 


[jeu]  =  jr 


{E}' 


[A41000KR}  +  V3  [A41 100]{R} 


+  {R}'[A500001]{A} 


+  {P }T ( [A31 00 j  +  [A3001 ] ) 


[JEE]  =  {R}T([A3000]  +  V3 [A3010]) 


+  ^  (R}T 


[A40010](U}  +  v  3 [ A40 1 1 0 ] { U } 


[JEP]  =  f  -  {U}T 


[A3001 ]  +  [A3100] 


+  {U}  [A500001]{A} 
IP 


P 

Jj+1 


Jj+1 


[JPR] 

[JPU] 

[JPE] 


-(y-1){E}T[A3000]  +  {U}T [ A40000 3 { U } 


r  t  1  p 

[(y-1){R}  [A40000]{U}j 
[-  (y-l){R}T[A3000] 


2 
P 

j+1 
P 

j+1 


.  p 

-j+1 


[JPP]  =  [A200] 


Note  in  the  formation  of  equations  79-81,  that  the  Boolean  indices 

in  the  various  A-matrices  are  permuted  to  facilitate  differentiation  of 

each  expression  by  the  last  right  contraction  matrix.  Note  also  the 

considerable  commonality  pervading  the  Jacobian  construction.  The 

elements  of  each  Jacobian  are  formed  on  the  finite  element  domain  R1 
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using  the  element  matrices  listed  in  Appendix  A,  and  then  assembled 
into  the  global  form  using  the  operator  Se.  In  actual  practice,  the 
column  matrix  {5QI}  is  ordered  on  degrees  of  freedom  at  a  node,  e.g., 
{...,  <5Rj,  <SUj ,  6Ej,  6Pj,  <5Rj+i ,  Hence,  the  global  Jacobian  [J], 

equation  27,  is  block  tri-diagonal,  using  the  linear  (k  =  1)  finite 
element  formulation,  and  block  penta-di agonal  for  quadratic  (k  =  2) 
elements. 

For  the  conservative  form  of  the  governing  differential  equation 


system  7-10,  identify  the  volume-specific  dependent  variables  m  =  pu  and 
g  h  pe,  yielding 

l(p)  =  f£  [m  +  mA]  =  0  (82) 

L(m)  =  fjr  +  ^  [m2/p  +  p]  +  ir/pWax  =  0  (83) 

L(g)  =  ft  +  L(mg/p  +  mp/p)]  +  m/p(g+p)nA/3x  =  0  (84) 

L(p)  =  p  -  (y-1)  [g  -  m2/2p]  =  0  (85) 


For  equations  82-85,  A  is  defined  as  An  A(x).  The  omnipresence  of  in¬ 
verse  density  can  be  eliminated  by  introducing  the  (contravariant)  con¬ 
vection  velocity  (u)  v  =  m/p;  hence, 

L(p)  -  ft  +  h  +  PvAJ  =  0  (82A) 
L(m)  =  +  ^  [vm  +  P]  +  mv;>A/^x  (83A) 
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L(g)  =  ft  +  Jx  ^v9  +  vpl  +  v(?+p)9A/9x  =  0 


(84A) 


L(p)  =  P  -  (y-1)  [g  -  \  vm]  =  0  (85A) 

The  finite  element  solution  algorithm  statement,  equations  24-26, 
applied  to  equations  82A-85A,  yields 

{FR}  =  ([A200J  +  vi[A210]{R}j+1 

+  ^  |"[A210]{M}  +  v ^ { V }T C A30 1 1  ] { R } 

+  {A}T([A41000]{V}){R}1  =  {0}  (86) 

Jj+1  ,j 

{FM}  =  ([A200]  +  V2[A210]){M}:+1 

+  -mT[A3010]lM}-  [ A210] { P} 

+  v? { V  >T [ A301 1 ] (M }  +  {A}T([A41000]{V  >) {M } 

Jj+l,j 

=  (0)  (87) 

(FG>  =  ([A200]  +  v13[A210]){G}j:+1 

+  y  {V)T[A3010]{G}  +  {V)T[A3010]{P} 

+  v  ?  {  V  }T  £  A30 1 1 3  { G }  +  {A}T([A41000]{V}){GP}1 

Jj+1  ,j 

=  (0)  (88) 
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{FP}  =  [A200] {P }  -  (y-1)  [A200] {G}  -  j  {V}T[A3000]{M}  =  {0}  (89) 


Note  that  equations  86-89  are  rather  less  complicated  in  appearance  than 
equations  70,  75-77,  and  in  conformation  with  the  multi-dimensional 
generalized  coordinates  formulation,  see  equation  41.  The  definition 
of  convection  velocity  has  been  utilized  for  the  dissipation  (v2)  term 
in  equation  86  to  enhance  overall  uniformity,  and  {GP}  =  (G  +  P}. 

The  construction  of  the  Jacobian  contributions  for  the  Newton  iter¬ 
ation  algorithm,  equations  27-29,  proceeds  directly.  Recalling  that 
v  =  m/p,  and  using  the  chain  rule  as  required. 


3{FR> 

w 


=  [JRR]  =  [A200]  +  vJ[A210] 


[JRM]  =  ^  [[A210]  +  vJ[A211]  +  (A}T[A3100] 


j+1 


[JRG]  =  0 


[JRP]  =  0 


(90) 


w  s  «*] 


At  1 

m 

CSJ 

[p2J 

{M)T[-  [A3010]  +  v|[A31 10] 


-P 

-j+1 


+  [A40001 ]{A} 

[JMM]  =  [A200]  +  v2[A210] 

+  y  {V}T[[A3010]  +  v2[A3011]  +  C A40001  ] { A}J ^ 

+  M_Lj{M}T|[A30l0]  +  v2[A3110]  +  [A40001  ]  { A}J  ^ 


lP 
j+1 
P 

j+1 


[JMG]  =  [0] 


[JMP]  =  [A210] 
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(91) 


=  [JGR]  =  -  ^(JyJ  (G}T([A3010]  +  Vj[A31 10]  +  [A40001]{A}) 


[JGM]  =  y 


1 

(.Pi  L 


T  "|P 

+  {P}  [A3010] 

Jj+1 

{G}T ( [A301 0]  +  v![A31 10]  +  [A40001 ]{A}) 


+  {P}  [A3010] 


j+1 


[JGG]  =  [A200]  +  VJCA210] 


+  y{V}T 


[A3010]  +  V3[A3011]  +  [A40001 ]{A) 


[JGP]  =  f{V}T 


[A3010]  +  [A40001]{A> 


P 

Jj+1 


-|P 

j+1 


(92) 


3{FP) 

9{R> 


=  [JPR] 


y-1 


{M  } 1  [A3000] 


-P 

-j+1 


[J  PM]  =  ^y 


|^{  V  }  [A3000]  + 
[JPG]  =  - (y-1 )[A200] 


r4-l  (M}T[A3000] 

.  p  ; 


P 

Jj+l 


[JPP]  =  [A200] 


(93) 


In  equations  91-93,  the  superscript  bar  on  m  and  p  indicate  the  assembly 

of  elemental  averages  of  the  variable,  i.e.,  me  =  {A10}^(M}e,  Pp  = 
t  e 

{A1 0}  {R)e.  The  various  defined  matrices  are  listed  in  Appendix  A, 

for  k  =  1  and  2  in  equation  23. 


2.  Numerical  Results,  Riemann  Shock  Tube 

The  primary  requirement  for  the  computational  tests  is  to 
evaluate  the  various  dissipation  constants  within  the  finite  element 
algorithm  for  compressible  flowfields  possessing  sharp  gradients  and 
shocks.  An  equally  important  requirement  is  to  assess  construction  of 
the  Jacobian;  specifically,  is  the  favorable  quadratic  convergence  rate 
for  the  Newton  iteration  algorithm  retained  in  the  present  formulation? 
Two  test  cases  meeting  these  requirements  are  the  Riemann  problem  and 
off-design  (shocked)  flow  in  a  nozzle.  Following  some  numerical  exper¬ 
imentation,  the  nonconservative  formulation,  equations  70,  75-77,  was 
relegated  to  history  in  favor  of  the  conservative  algorithm  form. 

The  Riemann  problem  corresponds  to  the  one-dimensional,  ex¬ 
ploding  diaphragm  shock  tube  experiment  (see  Shapiro,  ref.  13.,  page 
1007),  with  uniform  sound  speed  in  both  the  high  and  low  pressure 
chambers.  Upon  (computational)  removal  of  the  diaphragm,  a  shock  wave 
propagates  into  the  initially  low-pressure  region,  and  a  rarefaction 
wave  propagates  in  the  opposite  direction.  The  more  interesting  case, 
corresponding  to  different  initial  sound  speeds  in  the  two  chambers, 
has  been  recently  examined  in  considerable  detail  by  several  researchers 
in  computational  methodology.  Particular  reference  has  been  made  to 
establishment  of  the  "ultimate"  flux-corrected  conservative  difference 
scheme,  cf . ,  Sod  (ref.  14),  Van  Leer  (ref.  15),  Zalesak  (ref.  16). 

These  methods  generally  employ  a  Lagrangian,  or  mixed  Lagrangian- 
Eulerian  formulation,  in  distinction  to  the  implicit,  strictly  Eulerian 
finite  element  algorithm,  which  is  not  at  all  "hard-wired"  for  the 
problem. 

The  selected  Riemann  problem  specification  corresponds  to  that 
of  Sod  (ref.  14),  wherein  a  unit  duct  length  is  uniformly  subdivided  by 
100  nodes.  The  diaphragm  is  located  at  node  50,  with  the  initially 
high  pressure  region  occupying  the  left  half.  For  the  finite  element 
simulation,  the  linear  (k  =  1)  algorithm  employed  99  uniform  elements, 
while  for  the  quadratic  (k  =  2),  the  last  node  was  deleted  and  the 
discretization  contained  49  elements.  The  initial  conditions  are 
u  (x)  =  0,  p  =  p  =  1  for  x  <  0.5,  p  =0.1,  and  p  =  0.125  for  x  >  0.5, 


and  y  =  1.4.  Figure  6  shows  the  base  case  k  =  1  algorithm  solution,  at 
t  =  0.14154s,  for  v2  =  1//15  and  v*  =  0,  where  ot  denotes  dependent  vari¬ 
able  (number)  and  the  dissipation  parameter  indices  (1,2)  are  now  super¬ 
scripts.  The  shock  is  located  at  node  75,  and  the  rarefaction  wave  is 
centered  on  node  62.  The  symbols  correspond  to  nodal  values  of  (Ql)  ,  while 
the  various  solid  lines  bound  the  inviscid  solution  (ref.  15)  for  the  shock, 
rarefaction  wave  and  the  interspersed  plateaus.  Considering  the  simpli¬ 
city  of  the  k  -  1  algorithm,  these  results  are  quite  accurate,  except  for 
the  modestly  sloped  plateaus  in  density  and  momentum,  and  the  overshoot 
for  all  (Ql)  on  the  high  pressure  side  of  the  shock. 

Following  some  experimentation,  the  improved  results  shown  in  Fig.  7 
were  obtained  for  v*  =  1//15  (3/8,  0,  1/4)  and  v2  =  1//15  (3/4,  2,  1), 
for  1  <  a  <  3.  The  plateau  in  momentum  and  density  behind  the  shock  is 
nearly  planar,  no  overshoot  occurs  in  any  (Ql),  and  the  shock  is  inter¬ 
polated  across  approximately  two  elements.  The  second  density  plateau 
is  also  modestly  improved,  but  the  "trashiness"  in  {Ql }  on  x  <  0.4  has 
been  somewhat  aggravated.  Since  the  velocity  is  zero  in  this  region, 
the  algorithm  is  operating  without  any  dissipation,  and  v1  >  0  in  the 
absence  of  dissipation  is  well  known  to  induce  this  type  of  leading  phase 
error. 

For  comparison.  Figure  8  is  a  reproduction  of  the  results  of  Van 
Leer  (ref.  15,  Figure  6),  for  this  problem  specification,  as  obtained 
using  the  MUSCL  algorithm.  The  upstream  density  plateau  of  this  algo¬ 
rithm  is  better,  but  the  results  in  the  shock  vicinity  are  nominally 
comparable.  The  solution  parameters  of  velocity  and  internal  energy 
are  rather  better  descriptions  for  algorithm  performance,  see  the  two 
right  side  plots  in  Figure  8.  Figure  9a  shows  these  data  computed  from 
the  k  =  1  solution  of  Figure  7.  Comparing  to  Figure  8,  the  shock  sharp¬ 
ness  is  comparable,  as  is  the  high  temperature  plateau  behind  the  shock. 

The  low  temperature  plateau  behind  the  rarefaction  wave  is  modestly 
less  planar.  For  comparison.  Figure  9b  shows  these  identical  data  as 
obtained  using  the  Crank-Nicolson  implicit  finite  difference  equivalent 
of  the  k  =  1  algorithm,  with  v2  =  1//T5  and  v1  =  0  (by  definition). 

This  ad  hoc  modification  is  again  verified  to  degrade  solution  accuracy. 
Upstream  overshoot  is  excessive,  the  shock  is  interpolated  across  six 


Figure  6.  Linear  Finite  Element  Algorithm  Solution,  Shock  Tube  Problem, 
t  =  0.14154,  v1  =  0.,  v2  =•  1//15. 


Figure  8.  Solution  for  Shock  Tube  Problem  Generated  By  the  MUSCL  Code 
Reported  by  Van  Leer  (ref.  15),  Courant  No.  =  0.9,  Ax  =  0.01 
t  =  0.14154. 


a)  Linear  Finite  Element,  Solution  of  Figure  7. 


b)  Finite  Difference,  =  l//T5\ 

'  a 

Figure  9.  Finite  Element  and  Finite  Difference  Algorithm  Solution  Com 
parisons.  Shock  Tube  Problem,  t  =  0.14154 
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difference  cells,  and  the  high  temperature  plateau  is  barely  recognizable 
and  underpredicted. 

The  performance  of  the  quadratic  (k  =  2)  algorithm  formulation 
is  an  improvement,  in  accord  with  previous  numerical  experience.  Figure 
10  shows  the  k  =  2  solution  field  {Q I }  at  t  =  0.14154  s,  using  the  linear 
analysis  dissipation  levels  =  1//1"5,  v*  =  0.  The  solution  is  devoid 
of  the  plateaus,  and  the  shock  is  smeared  over  several  elements,  in¬ 
dicative  of  excess  diffusion.  Figure  11  illustrates  the  improvements 
accrued  to  use  of  v*  =  1//T5  {1/4,  3/4,  1/2}  and  v*  =  0.  The  shock  exists 
across  only  one  element  and  the  density  and  momentum  plateaus  are  clearly 
evident.  A  modest  level  of  solution  "trashiness"  is  evidenced  in  both 
extreme  solution  regions  wherein  the  dissipation  level  is  zero  (due  to 
zero  convection  velocity).  Figure  12a  shows  the  velocity  and  internal 
energy  fields  computed  from  this  quadratic  solution,  and  the  fidelity  is 
generally  excellent.  In  this  instance  (the  quadratic  element  algorithm), 
the  reduction  of  the  initial-value  matrix  to  the  Crank-Nicolson  finite 
difference  equivalent  does  not  adversely  affect  solution  accuracy,  Figure 
12b,  except  in  the  high  temperature  plateau  behind  the  shock. 

The  construction  of  equation  system  Jacobian,  as  presented  for 
both  k  =  1  and  k  =  2,  retains  the  favorable  nominally  quadratic  rate  of 
convergence.  Table  4  presents  the  extremum  elements  of  {<5QI } ,  for  a 
typical  integration  step  involving  three  iterations,  and  the  nodal 
location  of  these  extrema.  Convergence  is  at  least  quadratic,  including 
the  algebraic  pressure  equation.  The  sharply  defined  solutions,  as  ob¬ 
tained  using  minimal  levels  for  v2 ,  typically  required  250  iterations  to 
reach  t  =  0.14154  s.  This  corresponds  to  an  extremum  Courant  number 
(C  =  |u+a|At/Ax)  of  approximately  0.35.  No  attempt  was  made  to  extremize 
Courant  number  during  this  study. 

3.  Numerical  Results,  Shocked  Nozzle  Flow 

The  results  for  prediction  of  an  off-design  nozzle  problem  confirm 
these  data  for  dissipation  parameter  levels.  The  test  case  initial  con¬ 
dition  corresponds  to  a  subsonic-supersonic  expansion  on  0.75  s  M  s  1.35, 
followed  by  a  normal  shock  and  subsonic  expansion  on  0.75  <  M  <  0.7.  The 
values  of  .  ,  .,u  and  ,>e  were  specified  at  the  inlet,  and  p  was  computed  from 
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Figure  11.  Quadratic  Finite  Element  Algorithm  Solution,  Shock  Tube 
Problem,  t  =  0.14154,  v*  =  0.,  v2  =  1 //TF{ 1/4,  3/4,  1/2}. 
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TABLE  4 


CONVERGENCE  IN  {SQI }  FOR  RIEMANN 
SHOCK  TUBE  SIMULATION 


Iteration 

HBMlabtiiiE  - 

Node 

{<$M}max 

Node 

1 

-0.87325E-02 

70 

-0.20780E-01 

70 

2 

0. 13072E-02 

69 

0.31246E-02 

69 

3 

0.21915E-03 

69 

0.69326E-03 

69 

{SGlmax 

Node 

Node 

1 

-0.23351E-01 

70 

-0.75154E-01 

70 

2 

0.34658E-02 

70 

-0.51826E-02 

71 

3 

0.46731E-03 

70 

-0.50114E-03 

68 

the  equation  of  state.  The  subsonic  outlet  conditions,  applied  at  the 
end  of  a  uniform  cross-section  extension  of  the  nozzle,  were  specified 
pressure  and  vanishing  normal  derivative  for  p,  pu  and  pe.  For  the  test, 
the  exit  pressure  was  raised  by  15%,  such  that  the  shock  must  move  up¬ 
stream  into  a  previously  supersonic  region  of  flow.  At  all  points,  up¬ 
stream  of  the  translating  shock  location,  the  supersonic  flow  must  remain 
undisturbed  from  the  initial  conditions.  Figure  13  shows  the  k  =  1 
algorithm  solution  for  {QI }  and  Mach  number  distribution  at  t  =  2.4s, 
compared  to  the  initial  conditions  shown  as  a  solid  line.  The  vari¬ 
ables  upstream  of  the  new  shock  location  are  unaltered,  and  the  downstream 
distributions  are  smooth.  The  specified  levels  of  dissipation  parameters, 
set  at  one-half  those  of  the  Riemann  problem,  were  adequate  to  suppress 
overshoot  and  to  produce  a  sharply  defined  shock.  Figure  14  presents  the 
similar  comparison  for  the  quadratic  algorithm  solution,  as  obtained 
using  one  half  the  dissipation  parameter  levels  of  the  Riemann  specifi¬ 
cation.  The  results  are  essentially  indistinguishable  from  those  of 
Figure  13. 


► 


48 


X  COORDINATE  x  COORDINATE 


Figure  13.  Computed  Solution  for  Shocked  Flow  in  a  Variable  Area 

Duct,  Linear  Finite  Element  Algorithm, - Initial  Condition. 
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Figure  14.  Computed  Solution  for  Shocked  Flow  in  a  Variable  Area 

Duct,  Quadratic  Finite  Element  Algorithm, - Initial  Condition. 
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Two-Dimensional  Flow  Formulation 


The  extension  of  the  formulation  of  the  dissipative  finite  element 
algorithm  to  a  two-dimensional  flow  description  in  generalized  coordi¬ 
nates  is  direct.  Identify  the  volume-specific  dependent  variable  set  p, 
m.  =  pu^ ,  1  <  i  <  2,  and  g  =  pe,  and  the  parameters  Oj j  and  q j.  Recall 
the  definition  of  the  convection  velocity,  as  expressed  in  contravariant 
scalar  components  equation  38.  For  a  multi-dimensional  problem,  the 
dissipation  parameter  B'  is  a  vector,  with  contravariant  scalar  compo¬ 
nents  v^Ag,  see  equation  24,  where  a  denotes  initial-or  boundary-value 
definition.  Finally,  as  before  for  compressible  flow,  B3  =  0. 

Equation  26  expresses  the  non-linear  algebraic  equation  system  resul¬ 
ting  from  substitution  of  equations  1-7  into  equation  24.  Denote  the  dis¬ 
crete  dependent  variable  set  {Q}^  =  {R,  M(I),  G,  P,  S(I,0),  Q(I)).  The 
respective  algorithmic  statements  {FI}  in  the  generalized  coordinate 
description  become,  see  equations  37-41, 

{FR}  =  ({DET}T[B3000]  +  vj j{ETAK0}T ( [B40K00] fDET}) ){R} j+] 

+  A*  -{ETAKI}T[B30K0]{Mn  +  v*j{ETAKJ}T[B30KL]{MBARL}  (94) 

L  j+l,j 

{ FMI }  =  ( {DET }T[B3000]  +  vjj{ETAKJ}T([B40K00](DET}) ){MI}j+1 

+  ^  -{UBARK}T[B30K0]{MI}  -  {ETAKI}T[B30K0]{P} 

+  {ETAKJ}T[B30K0]'SIGIJ} 

+  Vj  .{ETAKJ }T ( [B40KLG] {MI } ) { UBARL } 

+  v?ij{ETAKJ}T([B40KL0]{UBARL}){MI}  .+1  j  (95) 
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{FG}  =  ({DET}T[B3000]  +  vJj{ETAKJ}T[B40K00]{DET})){6}j+1 

+  ^[-{UBARK}T[B30K0]({G}  +  IP})  +  {ETAKL}T[B30K0]{QBARL} 

+  {ETAKJ}T([B40K00]{UL}){SIGLJ} 

+  vfj  iETAKJ  }T  ( [  B40  KLO]  { G} )  {  UBARL } 

+  v2J{ETAKJ}T([B40KL0]{UBARL}){G}  J+1  1  (96) 

{ FP }  =  {DET }T[B3000]{P}  -  (y-1 ) ({DET}T[B3000]{G} 

-  i{DET}T([B40000]{UJ}){MJ})j+1  (97) 

A  few  comments  on  notation  in  equations  94-97  are  appropriate. 

The  matrix  B-prefix  denotes  the  expressions  obtained  by  integrals  over  a 
two-dimensional  elemental  domain,  i.e.,  R*.  The  defined  element  matrices 
are  listed  in  Appendix  B  for  k=l ,  equation  23.  The  concept  of  the  Boolean 
indices  (0,1)  is  generalized  to  replace  1  by  K,  indicative  of  the  corre¬ 
sponding  scalar  component  of  V>  in  the  coordinate  system,  hence, 

1  <  K  <  2.  In  all  terms,  the  discrete  indices  J,  K,  L  occurring  in  both 
matrix  and  variable  (FORTRAN)  names,  are  tensor  summation  indices  with 
range  1  <  (J,  K,  L)  <  2.  The  index  I  is  a  free  index  denoting  Cartesian 
scalar  components  of  m^ ,  i.e.,  {MI}.  As  noted,  the  expansion  coefficient 
.  equation  24,  is  expressed  in  Cartesian  scalar  components  Vj,  with 
distinct  values  for  each  dependent  variable  and  each  term.  The  arrays 
{DET}  and  {ETAKJ}  contain  nodal  values  of,  respectively,  the  determinant 
J  of  the  forward  transformation  *.  =  x. (n 

inverse  transformation,  see  equations  33-31.  The  nodal  values  of  the 
contravariant  scalar  components  of  the  solution  vectors  u^,  and  m^  and 
qR  are  denoted  {UBARK},  {MBARK}  and  {QBARK}.  The  elements  of  {SIGIJ} 
are  nodal  values  of  the  stress  tensor.  The  corresponding  algorithm 
statement,  using  the  homogeneous  form  of  equation  4,  is 


j)  and  elements  9r^/3xj  of  the 


{ FS I J }  =  {DET}T[B3000]{SIGIJ} 


-  y 


{ETAKJ}  ([B400K0]{DET }){UI} 


+  {ETAKI }  ( [ B400 KO ] { DET } ) { U J } 


-  -^<5  j  j  {  ET AKL  }T  ( [  B400K0]  { DET } )  {  UL } 


j+1 


(98) 


where  6^  is  the  Kronecker  delta.  The  corresponding  expression  for 
the  definition  of  energy  flux  q.,  equations  5-6,  is, 

J 


(FQI)  =  {DET}T[B3000]{QI } 


-  k  {ETAKI }T([B40K00]{DET}) ({GSR}  -  {UIUI}) .+1  (99) 


In  equation  99,  the  elements  of  (GSR)  are  nodal  values  of  e  h  g/p, 
while  those  of  {UIUI}  are  nodal  values  of  specific  kinetic  energy  iu.u^. . 

In  all  equations,  the  notation  {  •  }j+1  denotes  {  •  }P+^  -  {  •  }j,  and 
]j+l  j  denotes  evaluation  of  the  argument  at  tj+-|  and  tj  followed  by  addi¬ 
tion. 


As  can  be  observed,  the  finite  element  two-dimensional  Navier-Stokes 
algorithm  statement  in  generalized  coordinates,  equations  94-99,  is  con¬ 
siderably  more  involved  than  the  predecessor  one-dimensional  algorithm. 

By  the  same  token, equations  94-99  also  completely  express  the  three- 
dimensional  algorithm  by  the  exchange  of  B-prefix  matrices  with  C  and  ex¬ 
tension  of  the  discrete  tensor  index  range  to  1  <  (J,K,L)  <  3. 

The  remaining  step  in  the  algorithm  formulation  is  construction  of 
the  tensor  matrix  product  form  of  the  Jacobian, 


[J({FI}]  = 


9{FI } 

aTqoT 


(100) 
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equation  29.  The  construction  and  subsequent  solution  procedure  was 
outlined  in  Section  II. 3.  The  lead  term  in  each  of  the  equations  is 

[DET]T[B3000] {QJ } j+] 

Hence,  the  corresponding  component  of  [J]  is 

[J]i  =  Sr  =  ^ET}T[B3000]6I0  (101) 

For  illustration,  assume  an  affine  rectangular  Cartesian  coordinate 
transformation,  det[J]  =  1,  hence  {DET }"*"  =  {ONE}^,  the  array  of  unit 
values.  By  direct  substitution, 

(0NE}T[B3000]  =  [B200]  (102) 

and  the  matrix  whose  tensor  product  yields  [B200]  is  simply  [A200],  see 
equations  49-51.  Recalling  that  global  expressions  are  obtained  by  the 
assembly  operation  Sg  over  element-level  operations,  equation  101  is  for¬ 
mulated  as 

{DET }T[B3000]  =  Se  Ae{DET}^[B3000]e  (103) 

=  Se  Ai{DET)T[A3000]e&A2{DET}T[A3000]e 

Recall  the  superscripts  on  Agl  denote  the  appropriate  one-dimensional  ele¬ 
ment  measure,  e.g.,  length  and  width  of  a  rectangular  element. 

The  second  term  in  the  Jacobian  for  p  and  a  is 

[J] 2  =  vjj{ETAKJ}T[B30K0] 

where  I  equals  1  or  4  for  two-dimensional  flow,  and  summation  is  implied 
for  repeated  indices.  By  the  direct  extension  of  operations  yielding 


equation  103, 


vi 1{ETAKJ}T[B30K0]  =  SQ  A  uJ.{ETAKJ}I[B30K0]a 

1J  6[  c  Id  6  6 

=  Se  AevJj{ETAlJ}^[A3010]e 

®  AJvi j {EtA2J }TL A30 1 0 ]el  (104) 

As  before,  the  superscripts  on  A^  denote  coordinate  direction,  with  the 
corresponding  indication  in  the  elemental  arrays  {ETAKJ}.  The  element 
derivative  matrix  [A3010]  is  independent  of  a  and  identical  with  the  one¬ 
dimensional  form  given  in  Appendix  A. 

The  construction  of  certain  remaining  terms  in  the  Jacobian  involves 
differentiation  with  respect  to  the  parameter  u^  =  m^/o.  Using  the  chain 
rule  and  equation  38-39,  then 

U _  .  2 _  +  .3 _ <HMK)  +  2 _ 9 { Ok) 

TImIT  7m  :  tMK>  9  (M I  >  9{  UK}  9 (M I } 

=  ki +  de1 

JL  -  1  . 

tTRT  '  9TRT  $ 

k 

which  is  the  generalization  of  the  one-dimensional  construction.  Recalling 

again  that  all  operations  are  performed  on  the  elemental  level,  denoting 

D„.  as  the  element  average  of  det  J  ,  and  letting  K  signify  the  dis- 

KI  L  G 

Crete  free  index  corresponding  to  ,  the  non-empty  tensor  product 

Jacobians  for  equations  94-99  are: 

[JRR  ]  =  1DET  T[A3000]  +  v’lETAKJ^CASOlO] 


[JRMI]=-y  -{ETAKI  }T[A3010] 

+  vfj8KI(ETAig)T[A3011]J  (106) 

[JMIR]  =  [  (MI  }T[A3010] 

+  VIJ  {ETAKJ}T([A40110]  +  [A40101]){MI > 

[JMIMI]  =  {DET}T[A3000]  +  vJj{ETAKJ}T[A40100]{DET} 

+  -{UBARK}T[A3010]  -  Dki(MI}T[A3010]^ 

m 

T  [  (  )  ( 

+  Vjj  {ETAKJ }  [A40101]  +  [A40110]  { UBARK} 

*  0Hl| 

; 

[JMIKJ]  =  x(^Kj)  i  -{MJ}T[A3010] 

+  VjL{ETAKL}T(([A40101]  +  [A40110]){MJ}) 

[OMIP]  =  -y-{ETAKI)T[A3010] 

[JMISIJ]  =  ^[{ETAKJ}T([A40100]{DET})  (107) 

[JGR]  =  ^^)[(G+P}[A3010] 

+  {ETAKJ  }T([A40110]  +  [A40101  ])  {G} 
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[OGMI]  =x(DKl](lj  -^G+pJTCA30103 


+  v|l{ETAKL>T ( ([A40101]  +  [A40110]) {MI}) 

[JGG]  =  (DET }T[A3000]  +  vjj {ETAKJ }T[A3010] 

+  -  { U  BAR  K } T [ A3  010] 

+  v|j {ETAKJ }T( [A401 1 0]  +  [A40101 ] ) {UBARK}J 
[JGP]  =  -  4r(UBARK}T[A3010] 

[JGSI J]  =~  {ETAKJ}1  ( ( {UI }T[A500100] ) {DET}) 

[JGQI]  =  ^feLIj{ETAKL}T[A3010]  (108) 

[JPR]  =  -  2r  {DET}T([A40000]{MK} 

[JPMK]  =  ^  {DET}T([A40000]{UK})  +  fij  {DET }T  ( [A40000] {MK} ) 

[JPG]  =  -(y-l){DET}T[A3000] 

[JPP]  =  { DET }T [ A3000  3  (109) 

The  Jacobi  arts  for  {FSIJ}  and  {FQI }  are  constructed  in  the  manner  indenti- 
cal  to  equation  109. 

Observing  equations  94-109,  the  multi-dimensional ,  generalized  coor¬ 
dinate  formulation  is  considerably  more  detailed  than  the  one-dimensional 
form.  However,  the  calculus  and  algebra  procedures  build  directly  upon 
the  elementary  concepts. 


5 •  Numerical  Results 

The  numerical  evaluation  of  the  two-dimensional  algorithm  is  in  its 
beginning  stage,  and  only  cursory  results  are  completed.  A  critically 
important  confirmation  of  the  tensor  product  Jacobian  formulation  is  com¬ 
pleted  for  the  scalar  dependent  variable  and  the  rotating  cone  test  case. 
Figure  15  summarizes  the  significant  portion  of  the  k  -  1  algorithm  solu¬ 
tions  fields  at  the  quarter  turn.  The  results  in  Figure  15a)  were  pro¬ 
duced  using  the  time  split  formulation  (ref.  5  ),  which  may  be  considered 
as  exact.  Figure  15b)  shows  the  tensor  product  Jacobian  algorithm  results, 
generated  by  the  algorithm  specification  in  equations  42-45  for  a  two- 
dimensional  problem.  Agreement  is  within  the  +1  band  allowed  by  trunca¬ 
tion  of  floating  point  to  the  integer  data  used  for  output.  For  compari- 
son,Figure  15c)  is  the  solution  produced  by  the  tensor  product  algorithm, 
applied  to  equations  26-27  rewritten  to  account  specifically  for  linearity, 

[Ju]{Ql)j+1  »  (RHS) j  (110) 

where  { RHS> ^  is  the  time  level  j  evaluation  of  the  difference  in  the  ini¬ 
tial  value  and  convection  contributions,  see  reference  5,  equation  36. 

The  tensor  product  Jacobian  form  of  equation  110  is, 

[Ji]{Pl)j+1  =  (RHS)j 

[J2]{Ql}j+1  =  (P1}J+1  (111) 

where  (Pll  is  an  intermediate  solution.  Comparing  Figures  15a)  -  15c) 
confirms  this  procedure  introduces  significant  error  into  the  solution, 
in  distinction  to  the  accurate  representation  provided  by  the  (Newton) 
iterative  formulation. 

The  construction  of  the  generalized  coordinates  formulation  has  been 
validated  using  for  comparison  the  potential  solution  to  incompressible 
flow  about  a  cylinder.  For  this  case,  the  contravariant  components  of 
convection  velocity  correspond  (to  within  det[J]  =  r)  to  the  radial  and 
azimuthal  velocity  components.  The  construction  of  { FR)  and  (FMI),  and 
the  corresponding  Jacobian  tensor  product  matrices  were  also  verified. 
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Figure  15.  Confirmation  of  Tensor  Matrix  Product  Jacobian  Formulation, 
Linear  Finite  Element  Algorithm,  Rotating  Cone  Test  Case, 
One-Quarter  Turn. 
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6.  Continuity  Constraint  Formulation 


An  important  aerodynamics  problem  corresponds  to  low  subsonic 
(M  <  0.3)  flows,  wherein  the  fluid  is  essentially  constant  density.  As 
a  consequence,  the  time  derivative  term  is  lost  in  the  continuity  equa¬ 
tion  1.  Correspondingly,  in  many  of  these  aerodynamic  applications,  the 
steady  flow  is  dominantly  unidirectional,  amenable  to  prediction  using  a 
viscous  marching  procedure,  typically  termed  parabolic  Navier  -  Stokes. 
Equations  11-18  provide  the  appropriate  differential  equation  system. 

The  basic  requirement  is  to  develop  and  evaluate  a  numerical  algorithm 
for  enforcement  of  the  continuity  equation  upon  solution  of  the  parabolic 
partial  differential  equation  system.  Specifically,  under  the  order  of 
magnitude  analysis  applied  for  generation  of  parabolic  Navier  -  Stokes,  the 
non-parabolic  continuity  equation  governs  first-order  effects,  while  the 
parabolic  transverse  momentum  equations  describe  second  and  higher-order 
effects. 

The  theoretical  concept,  borrowed  from  the  variational  calculus,  is  en¬ 
forcement  of  the  continuity  equation  (solution)  as  a  differential  constraint 
on  solution  of  the  transverse  momentum  equations.  This  measure  for  continuity 
must  be  global,  i.e.,  span  Rn,  and  the  deviation  will  vanish  as  continuity 
becomes  satisfied.  An  appropriate  global  measure  of  the  deviation  from 
the  solution  of  equation  11  is  the  harmonic  function  <J>,  satisfying 

LW  -°  (m) 

J  J 

subject  to  the  boundary  conditions  on  9R, 

*U)  =  =  0  (113) 

and  setting  <t>  =  0  at  least  at  one  location  on  9R.  As  equation  11  be¬ 
comes  satisfied  as  a  differential  constraint,  equations  112-113  become 
homogeneous  and  the  solution  vanishes.  Hence,  L(p^)  in  equation  24 
is  replaced  by  L(<t>h),  the  discretized  statement  of  equations  112-113. 

Assessment  of  algorithm  performance  focuses  on  determination  of  the 
expansion  coefficient  p>3,  as  well  as  implementation  of  the  algorithm  se- 
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quencing.  Following  consequential  numerical  experimentation,  f?3  =  At, 
the  downstream  marching  step,  was  determined  suitable,  coupled  with  a 
procedure  for  summing  sequential  contributions  to  the  constructed  deriva¬ 
tive  matrix.  Specifically,  considering  for  exposition  boundary  layer  flow, 
the  most  elementary  form  of  parabolic  Navier  -  Stokes,  the  algorithm  state¬ 
ment  (equation  26)  for  the  u2  momentum  equation  is 

{FU2}  e  {U1 }T[A3000]({U2}j+i  -  (U2> 

+  At({GU2}P+1  +  (GU2}J)  =  {0}  (114) 

where  {GU2}  contains  the  contributions  due  to  spatial  distributions  of 
transverse  convection,  pressure  and  viscosity,  as  well  as  the  continuity 

X  L 

constraint.  At  the  pin  iteration  for  step  j+1,  following  solution  of 
L(4>*1)  for  {<j>}j+.|,  the  contribution  to  the  continuity  constraint  is  evalu¬ 
ated  as, 

(GU20}P+1  =  [A2TO]{0) j+i  (115) 

This  contribution  is  accumulated  into  the  previous  evaluations,  yielding 

(GU21j*i  s!£,  {gu2*>!m  <"6) 

as  the  action  of  the  continuity  constraint  in  equation  114.  Hence,  each 
successive  determination  of  {0} corrects  the  action  of  all  previous  itera¬ 
tions,  such  that  -Ke},  where  |e|  >  0  is  an  acceptable  discrete  level 

of  computed  zero,  in  the  limit  as  p  increases  without  bound.  This  proce¬ 
dure  thus  admits  the  correct  (continuity  preserving)  solution  for  equations 
112-113,  i.e.,  0  *  0  everywhere.  Using  the  finite  element  algorithm  state¬ 
ment,  equation  24,  with  =  ft,  the  solution  statement  for  equations  112- 
113  is. 

{ FPH I }  =  [A2U]{0}j+i  +  ([A200]{U1}^+1  +([A201]{U2}P+1 )  (117) 

61 


r 


A  backwards  difference  formula  is  used  to  evaluate  the  elements  of 
{Ul}j+1.  The  various  A  prefix  matrices  are  familiar,  see  Appendix  A. 

7.  Numerical  Results 

A  definitive  test  case  for  the  continuity  constraint  algorithm  is 
provided  by  the  two-dimensional  boundary  layer  equations  for  laminar  flow 
in  zero  pressure  gradient.  The  freestream  level  of  u i ( xx )  remains  invari¬ 
ant,  yet  the  corresponding  plateau  in  u2(x2  ■>  f)  must  decay  uniformly  as 
x’J.  The  Poisson  formulation,  equations  112-113,  admits  the  required  glo¬ 
bal  coupling.  The  level  of  u2(x2)  ranges  over  five  digits  on  0  <  x2  ^  5, 
and  the  computational  solution  must  induce  a  vanishing  normal  derivative 
for  u2(x2  =  0),  see  equation  11,  since  9ui/3xi  =  0  at  x2  =  0. 

Using  a  M  =  32  non-uniform  discretization  of  R1,  the  results  generated 
using  the  conventional  finite  element  boundary  layer  algorithm  (ref.  3), 
and  the  continuity-constrained  parabolic  Navier-Stokes  algorithm  plot 
as  essentially  identical.  Table  5  summarizes  these  comparison  solutions, 
in  terms  of  coordinates  x2/iS  for  which  u^/u^  changes  by  an  order  of  mag¬ 
nitude.  The  agreement  is  excellent  over  the  entire  range,  as  is  the  ap¬ 
proximation  to  3u2/3x2  =  0  at  x2/6  =  0.  The  constraint  algorithm  converges 
to  { |  “QI  | }  <.  10’ 5  in  typically  4  to  5  iterations  per  step.  An  intrinsic 
error  measure  for  equation  112  is  the  energy  norm. 


E(4>,c|>) 


R1 


Ml. 

?x2  3xr 


(118) 


In  E(4>,4>)»  the  error  in  satisfaction  of  continuity  typically  decreases 
by  a  factor  or  two  for  each  iteration,  as  illustrated  in  Table  6  for  a 
typical  integration  step.  For  laminar  flow,  the  algorithm  maintains 

o  r 

E(4>,<j>)  <  10  for  iteration  convergence  to  t  *  10  . 

The  corresponding  levels  for  turbulent  boundary  layer  prediction  are 
-5  -4 

E (<t> ,d>)  <  10  for  c  ^  10  .  Figure  16  compares  the  continuity  constraint 
algorithm  results  to  data  (ref.  17)  for  the  Bradshaw  turbulent  relaxing 
flow  experiment.  This  solution  is  indistinguishable  from  the  direct 
boundary  layer  solution,  and  was  obtained  using  the  turbulence  kinetic 
energy-dissipation  function  closure  model  and  the  complete  equation  system 
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12-18.  The  continuity  constraint  algorithm  is  equally  applicable  to 
ducted  and/or  semi-bounded  flows,  and  should  find  considerable  use  in 
computational  aerodynamics. 


Table  5 

Transverse  Velocity  Distributions,  U2(x2)x103 
Laminar  Incompressible  Boundary  Layer 


Coordinate 

Boundary 

Continuity 

(x,/S) 

Layer  Solution 

Constraint  Solution 

0.0 

.0 

.0 

0.0009 

.0000011 

.0000001 

0.0021 

.0000054 

.0000030 

0.0035 

» 

.0000149 

.0000114 

* 

0.0095 

.00011 

.00011 

0.031 

.00116 

.00118 

0.10 

.0128 

.012/ 

0.67 

.139  * 

.138 

1.0 

.218 

.218 

Table  6 


Continuity  Constraint  Algorithm  Convergence 
Laminar  Incompressible  Boundary  Layer 


Iteration 

1 

2 

3 

4 


Energy  Norm,  E(^,4>) 
0.69432  E ( -9 ) 
0.30670  E { -9) 
0.15692  E(-9) 
0.08028  E ( -9 ) 


DOWNS1"  RE  A'-'  station.  XI /l 

Figure  16.  Computed  Boundary  Layer  Integral  Parameter  Distributions 

For  Bradshaw  Relaxing  Flow,  Linear  Finite  Element  Continuity 
Constraint  Algorithm,  Turbulence  Kinetic  Energy  Closure  Model. 
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SECTION  VI 


CONCLUSIONS  AND  RECOMMENDATIONS 

The  objective  of  this  research  project  is  to  derive  and  evaluate 
versatile,  accurate  and  efficient  numerical  al  g*' -i  thms  for  solution  of 
aerodynamic  flowfields  at  large  Reynolds  number.  Building  upon  work 
reported  in  reference  5,  the  concept  of  a  dissipative  finite  element  al¬ 
gorithm  has  been  refined  and  extended  to  solution  of  a  complete  equation 
set  in  aerodynamics.  The  test  case  numerical  results  are  highly  encour¬ 
aging  with  respect  to  shock  resolution  and  overall  performance,  utilizing 
both  a  linear  and  a  quadratic  finite  element  embodiment  of  the  theory. 

The  theoretical  formulational  statement  of  the  algorithm  has  been  extended 
to  a  multi -dimensional  description  in  generalized  coordinates.  The  key 
efficiency  feature  is  identification  of  the  tensor  matrix  resolution  of 
the  Jacobian  of  the  Newton  algorithm  for  this  statement.  Finally,  the 
concept  of  application  of  the  continuity  equation  solution,  as  a  differ¬ 
ential  constraint  on  the  momentum  equation  algorithm,  has  been  validated. 
The  formulation  is  directly  useful  for  viscous  marching  procedures,  and 
with  some  modifications  could  be  equally  useful  for  a  low  Mach  number 
Navier  -  Stokes  solution  algorithm. 

These  results  on  application  of  finite  element  solution  concepts  in 
computational  aerodynamics  are  highly  encouraging.  It  is  recommended 
that  the  multi -dimensional  computer  program  development  be  pursued  with 
vigor,  to  facilitate  performance  quantization  for  the  next  level  of  prob¬ 
lem  complexity. 
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APPENDIX  A 


FINITE  ELEMENT  ALGORITHM  HYPERMATRICES 
LINEAR  AND  QUADRATIC  BASIS  ON  ONE  DIMENSIONAL  SPACE 

1.  Linear  Finite  Element  Formulation 
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APPENDIX  B 


FINITE  ELEMENT  ALGORITHM  HYPERMATRICES 
LINEAR  BASIS  ON  TWO-DIMENSIONAL  SPACE 

Linear  Tensor  Product  Basis  (Ag  =  det  J): 
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